The present invention is directed to greatest common divisor algorithms and in particular to such an algorithm that uses multi-precision arithmetic and iteratively reduces the number of bits in each of the values it is testing either by performing an approximate division reduction step or by performing a novel inter-reduction step.
Modular arithmetic is used in many applications such as encryption algorithms and the generation of error correcting codes. When the modulus is a prime number, the operations of addition, subtraction, multiplication and division as well as the associative, commutative and distributive properties are defined over the set of numbers from zero to p. This set of numbers defines a Finite field modulo p, Fp.
GCD algorithms are commonly used to screen candidate large prime numbers. If a number passes the GCD test a more rigorous test is applied to determine if the number is prime. As set forth above, large prime number are important because they define Finite fields which are of interest for encryption purposes. One type of encryption algorithm encrypts data using exponentiation, relying on the inherent difficulty of the reverse of exponentiation, the discrete logarithm problem, to hold the data secure. Because encryption performed on a large Finite field is more secure than encryption performed on a small field, these algorithms work best when the exponentiation operations are carried out over a large Finite field. Thus, there is a need to identify large prime numbers. One problem with using large fields, however, is the size of the numbers being processed. Typical numbers used in data encryption have hundreds of bits. These numbers are too large to be easily handled by commonly available microprocessors that are limited to 32 or 64-bit arithmetic. As described below, calculations used in the encryption process are typically handled using multiprecision arithmetic.
FIG. 1 is a flow chart diagram of an algorithm that screens candidate prime numbers. The GCD algorithm is relatively quick and can find weed-out many candidate prime numbers having relatively small factors. The algorithm shown in FIG. 1 tests a number P to determine if it is prime relative to the first N thousand prime numbers where N is an integer. Prior to performing the algorithm shown in FIG. 1, an array THOUSAND_PRIMES [I] is prepared such that THOUSAND_PRIMES [1] contains the product of the first thousand prime numbers, THOUSAND_PRIMES [2] contains the product of the second thousand prime numbers and so on. The process begins at step 102 by obtaining the number P which is to be tested and the number N defining the number of thousand primes that P is to be tested against. At step 104 an index variable I is set to zero. At step 106, I is incremented by one and a value G is calculated where G is the greatest common divisor of P and THOUSAND_ PRIMES [I]. If, at step 108, G is greater than one then P is divisible by one of the factors of THOUSAND_PRIMES [I] and is not prime. Accordingly, control transfers to step 110 where the value FALSE is returned.
If, however, at step 108 G is equal to one then P is not divisible by any of the factors of THOUSAND_PRIMES [I] and control transfers to step 112. At step 112 the variable I is compared against N. If I is less than N control transfers to step 106 where I is incremented and the greatest common divisor of P and THOUSAND_PRIMES [I] is calculated for the incremented value of I. When, at step 112, I is equal to N, the number P has been tested against the N thousand prime numbers and has been found to be relatively prime to each of these products. Thus, at step 114, the algorithm shown in FIG. 1 returns a value TRUE indicating that P is a good candidate to be a prime number. After performing the test shown in FIG. 1, the candidate prime number may be processed using a probabilistic primality testing routine, such as the Miller-Rabin algorithm to determine if it is prime.
A key element of the program shown in FIG. 1 is the greatest common divisor (GCD) algorithm. FIG. 2 is a flow chart diagram which illustrates a greatest common divisor algorithm known as Euclid""s algorithm. The algorithm shown in FIG. 2 calculates the greatest common divisor of U and V where U is greater than V. The algorithm shown in FIG. 2 also calculates the inverse of V modulo U. This algorithm relies on the property that if U and V have a common divisor D so does U-V, U-2V and so on. Thus, using only subtraction, one can calculate the GCD of U and V. GCD algorithms run faster when they can combine multiple subtractions in a single step. When calculating the inverse, The GCD algorithm operates by maintaining the linear equations (1), (2) and (3)
U1P+U2X=U3xe2x80x83xe2x80x83(1)
V1P+V2X=V3xe2x80x83xe2x80x83(2)
T1P+T2X=T3xe2x80x83xe2x80x83(3)
Where Uxe2x89xa7V and U3 and V3 are initially assigned the values of U and V, respectively. Upon termination, U3=GCD(P, X)=1 and U2=Xxe2x88x921 MOD P. In general terms, GCD algorithms operate by repetitively reducing the number of bits in the larger value, U3, and switching the two values U3 and V3. Thus, the algorithm successively reduces the values of U3 and V3 while maintaining the linear equations. Because it maintains the values U2 and V2, the algorithm shown in FIG. 2 not only calculates the greatest common divisor of U and V but also calculates the inverse of V modulo U. That is to say, Vxe2x88x921 where V*Vxe2x88x921=1 modulo U. As described below, when U is a prime number, this inverse may be used to simplify division calculations performed within the Finite fields Fu. Furthermore, the variables U1, V1 and T1 do not need to be maintained because they can be determined from the other variables, for example, U1 can be determined from U2 and U3 by the identity U1=(U3xe2x88x92U2V)/U.
The Euclid GCD algorithm shown in FIG. 2 begins at step 210 by obtaining the values U and V. Next, at step 212 a temporary variable U3 is set equal to U and a temporary variable V3 is set equal to V. Also in step 212 variable U2 is set to zero and V2 is set one. Next, at step 214 the variable V3 is tested to determine if it is equal to zero. If not, step 218 is executed. This step calculates a value Q which is equal to the greatest integer less than U3/V3. Next, it determines a new value for V3 by storing the quantity, U3xe2x88x92Q*V3 into a temporary variable T3. Also at step 218, the process stores the value U2xe2x88x92Q*V2 into a temporary variable T2. This value is the new value for the variable V2. The value in T3 represents a reduced value of U3. Because V3 has not been reduced, however, T3 is less than V3. Consequently, the value in V3 is assigned to U3 and the value in T3 is assigned to V3. A similar operation is performed on the variables U2, V2 and T2. After step 218, control returns to step 214 where V3 is once again compared to determine if it is equal to zero. When V3 is equal to zero the process returns the values U3 and U2. The value U3 is the greatest common divisor of U and V while the value U2 is the inverse of V modulo U (i.e. Vxe2x88x921 MOD U).
As described above, a modular arithmetic operation that is used in many encryption algorithms is raising a large number M, representing a message, to a power E modulo P. In a typical encryption example the values M, E and P are each hundreds of bits in length. Operations on numbers of this sizexe2x80x94and especially division operationsxe2x80x94are difficult because common computer processing hardware supports at most sixty-four bit mathematical operations. Values greater than 64 bits in length are typically calculated using multiprecision arithmetic. A number, a, raised to an exponent e can always be calculated by multiplying that number by itself the number of time represented by the exponent, or in mathematical terms:
ae=a*a*a . . . e number of times.
Another method, which is significantly faster, is the multiply chain algorithm. In this case, let e=enxe2x88x921enxe2x88x922 . . . e1e0 be an n-bit exponent eixcex5{0, 1}, 0xe2x89xa6ixe2x89xa6nxe2x88x921 and enxe2x88x921=1. The algorithm starts with p1=a, then
Pi+1=Pi2 if enxe2x88x921-i=0 or a*pi2 if enxe2x88x921-i=1, where 1xe2x89xa6ixe2x89xa6nxe2x88x922.
Several methods are known in the art to reduce either the number of multiplications used to produce efficient exponentiation of the base value.
One method known to reduce the number of multiplication is the signed digit algorithm. Using this method, the exponent is represented as a string of bits comprising the values 0 and 1. Within the bit string, sequences (or xe2x80x9crunsxe2x80x9d) of 1""s are replaced by 0""s, with a 1 being placed in the next higher bit position to the most significant bit (MSB) position of the run, and xe2x80x9cxe2x88x921xe2x80x9d being inserted in the least significant bit (LSB) position of the run. By thus efficiently recoding the exponent bit string, the expected number of multiplication operations is reduced from n/2 to n/3.
Although the signed digit algorithm reduces the number of multiplication operations needed to calculated the exponential value Me modulo P, it replaces some of the multiplication operations by division operations. As described above, division may be implemented as the multiplication by the inverse.
FIG. 3 is a flow chart diagram which illustrates an exemplary method for calculating ME modulo P. The first step in this algorithm, step 310 gets the values P, M and E. P is a large prime number, M is the message to be encoded and E is an exponent value to be used in order to encode the message. The next step in the process is to determine the inverse of M modulo P (i.e. Mxe2x88x921 MOD P). This is accomplished by calculating the greatest common divisor of P and M. As P is a prime number by definition, the value U2 returned by the GCD algorithm is equal to Mxe2x88x921 MOD P. Next, at step 314, the signed digit algorithm is applied to the exponent E. Finally, at step 316, the multiply-square algorithm is applied to calculate ME modulo P. As described above, the multiply-square algorithm both multiplies and divides by M. The value Mxe2x88x921 MOD P is used to convert the division operation into a multiplication operation. This follows from the following equations (a/b*bxe2x88x921/bxe2x88x921) mod p=(a*bxe2x88x921)/(b*bxe2x88x921) mod p=(a*bxe2x88x921)/1 mod p.
The Euclid algorithm is not efficient for large numbers that require multiprecision arithmetic. Thus, there is a need for an efficient GCD algorithm that calculates inverses, and operates efficiently on numbers that require multiprecision arithmetic.
The present invention is embodied in a method for calculating greatest common divisors that uses an approximate division in its reduction step. The result of this approximate division is then compared to determine if the result is negative. If it is, then an alternate method is used to reduce the number of bits in the values while maintaining the linear formulas.
According to another aspect of the invention, if the first approximate division produces a negative result, the method applies a correction to the first approximate division to determine corrected values that have a reduced number of bits. If, during this correction step, the result is still negative, then yet another method is applied to reduce the number of bits in the values.
According to another aspect of the invention, the approximate division is applied only when the number of bits in the two values differ by at least a predetermined value. When the number of bits in the two values differ by less than this number, an alternative inter-reduction step is applied.